#Recuerden establecer el directorio de trabajo donde están sus datos
setwd("/Users/asernas/Desktop/2024/ASISTENCIAS/Genética\ de\ poblaciones\ 2024-1/LAB1/") Laboratorio 1 - Genética de Poblaciones
Introducción
La genética de poblaciones es una rama de la biología que estudia la variación genética de los individuos dentro de poblaciones. Parte de esta ciencia es obtener estimativas de la diversidad y la estructura genética de una o más poblaciones. Los datos o genotipos para el estudio de la genética poblacional se pueden obtener a través de marcadores moleculares. Los más utilizados recientemente, son los microsatélites o secuencias repetidas cortas (SSR) y los polimorfismos de un único nucleótido (SNP) (Jiménez, 2014).
Las principales estimativas utilizadas en el estudio de la genética poblacional determinan el nivel de polimorfismo, las frecuencias alélicas y las frecuencias genotípicas. El polimorfismo puede ser estimado a través del número de alelos por locus y por el porcentaje de loci polimórfico (Grünwald et al., 2015). A partir de las frecuencias alélicas y genotípicas, se pueden obtener otras estimativas como la heterocigosidad esperada y observada (Bedoya, 2012).
Para obtener estas estimativas existen diversos software y plataformas. Sin embargo, recientemente los autores recomiendan utilizar programas o lenguajes estadísticos más versátiles y de código abierto como R (Grünwald et al., 2015). El lenguaje estadístico R (R Core Team, 2013), ha permitido la creación de librerías o programas para el análisis de datos genéticos (Fuchs & Barrantes, 2015). Sin embargo, es necesario consultar o utilizar diversas librerías para realizar ciertos análisis. El objetivo de estos laboratorios será condensar la información de las librerías para trabajar con datos genéticos y presentar las herramientas para realizar análisis básicos de genética poblacional utilizando R.
Práctica 1. Lectura e importación de datos genéticos
Objetivo: Aprender a dar formato y a importar bases de datos de microsatélites y SNPs para realizar análisis genéticos utilizando distintas librerías.
R es un lenguaje y entorno para la computación estadística y confección de gráficos. RStudio es un entorno de desarrollo integrado elaborado para usar el lenguaje de programación R, el cual incluye una consola, editor de sintaxis que soporta la ejecución directa de código, así como herramientas para manejo de gráficos, historial y depuración. Ambos son de libre acceso y pueden ser descargados en sus sitios web oficiales. Para el uso de R con RStudio, es necesaria la descarga de tanto R (https://cloud.r-project.org/) como la versión gratis de RStudio (https://rstudio.com/products/rstudio/download/).
Una vez instalado tanto R como RStudio en su computadora, se procede a abrir RStudio y se crea un nuevo script de R, de la siguiente forma:
En ese script pueden digitar el código que se irá presentando en el manual, el cual se ejecuta presionando ctrl+Enter o ctrl+R. Una vez ejecutado el código, la salida del mismo aparecerá en la consola, la cual está ubicada en la esquina inferior izquierda de la ventana de R. En caso de que el código ejecutable dé como salida un gráfico, este aparecerá en “Plots”, en la esquina inferior derecha de la ventana.
En este manual, se colocará el código que deben ingresar en su script de RStudio. Este código usualmente es ejecutable, sin embargo, el texto después del símbolo numeral “#”, corresponde a comentarios y es texto no ejecutable. Estos comentarios se incluyen para agregar información, notas, instrucciones de código de uso personal para sus scripts.
1. Definición de directorio e instalación de paquetes:
Lo primero es definir el directorio de trabajo, que es la carpeta donde estarán los archivos para trabajar y donde se guardarán los archivos que se exporten de R. Para esto se utilizan los comandos:
Una vez definido, se procede a descargar, instalar y cargar los paquetes que se necesitan para trabajar, con los siguientes comandos:
## Ya se encuentran instalados en estos equipos, solo se deben cargar los paquetes
#install.packages( c("pegas", "adegenet", "poppr", "hierfstat","mmod", "magrittr", "ggplot2"), dependencies = TRUE)
library(pegas)
library(adegenet)
library(poppr)
library(hierfstat)
library(mmod)
library(magrittr)
library(ggplot2) 2. Librería pegas
Es un paquete que incluye funciones para leer, escribir, trazar, analizar y manipular datos alélicos y haplotípicos, además permite el análisis de secuencias de nucleótidos de población y microsatélites. Incluye análisis de coalescencia, desequilibrio de ligamiento, estructura de la población (Fst, Amova), HWE, redes de haplotipos, entre otros (Paradis, 2010).
Permite importar datos genéticos donde cada alelo de cada locus se almacena en columnas (número de columnas es igual al número de loci multiplicado por el nivel de ploidía). Para esto se utiliza la función alleles2loci(x, ploidy = 2, rownames = NULL, population = NULL) que transforma dichos datos en un objeto “loci”, un tipo de estructura de datos que utiliza pegas. Se le debe indicar los siguientes argumentos:
x |
Una matriz o un marco de datos donde cada columna es un alelo. |
ploidy |
Un número indicando el nivel de ploidía. |
rownames |
Un número entero que proporciona el número de columna que se usará para leer los nombres de los individuos. |
population |
Un número entero que proporciona el número de columna donde se indica el número de población |
Primero se debe importar el archivo .csv con los genotipos con el formato de un alelo por columna. Es importante revisar, y de ser necesario, editar los datos en caso de que exista algún punto “.” en el nombre de los SSR, es decir en el nombre de las columnas de los alelos. Es muy importante revisar el tipo de separación que tienen las columnas, ya sea por comas “,” o por punto y coma “;”, ya que Excel usa este último cuando se utiliza una coma como separador decimal.
Como ejemplo se utilizarán datos que corresponden a genotipos obtenidos con 20 SSR de muestras de maíz criollo de las regiones Brunca y Chorotega de Costa Rica. Dentro de cada región se definieron 10 sitios de muestreo en los que se colectaron mazorcas de distintos colores. La base de datos se llama maiz.csv
El siguiente código lee el archivo con formato csv y lo almacena como una tabla en el objeto data.maiz:
data.maiz <- read.csv2("maiz.csv")
head(data.maiz) id pop bnlg1046 bnlg1046.1 bnlg1191 bnlg1191.1 bnlg1194 bnlg1194.1
1 BE12.372 Brunca 173 173 197 199 188 188
2 BE12.381 Brunca 173 173 199 187 188 168
3 BE12.382 Brunca 173 173 189 189 168 168
4 BE12.383 Brunca 173 173 187 189 168 180
5 BE12.543 Brunca 196 196 187 219 186 186
6 BE12.544 Brunca 196 196 197 195 180 186
bnlg1265 bnlg1265.1 bnlg1449 bnlg1449.1 bnlg1523 bnlg1523.1 bnlg1740
1 196 221 87 148 191 191 123
2 215 233 87 148 NA NA 119
3 196 196 87 146 242 246 170
4 196 235 87 87 191 191 123
5 213 217 148 148 242 187 119
6 241 243 148 148 242 187 119
bnlg1740.1 bnlg1890 bnlg1890.1 bnlg1917 bnlg1917.1 bnlg2305 bnlg2305.1 phi015
1 123 161 204 NA NA 169 183 94
2 170 161 161 NA NA 183 183 94
3 170 94 94 98 98 169 183 97
4 119 94 94 98 98 169 169 94
5 119 94 125 98 128 189 189 97
6 119 94 94 98 98 NA NA 81
phi015.1 phi031 phi031.1 phi064 phi064.1 phi072 phi072.1 phi078 phi078.1
1 94 185 222 84 104 141 141 120 123
2 97 222 222 84 104 141 149 120 125
3 97 185 185 104 100 141 150 120 124
4 97 185 185 104 100 141 150 124 124
5 81 222 222 76 82 150 162 120 123
6 92 222 188 76 76 150 162 120 120
phi089 phi089.1 phi093 phi093.1 phi109188 phi109188.1 phi116 phi116.1
1 87 95 286 286 163 166 151 162
2 95 95 286 286 163 166 162 162
3 87 87 286 283 166 166 162 162
4 95 95 286 283 166 166 162 162
5 95 95 288 288 166 166 151 162
6 95 95 288 288 166 166 162 162
phi96100 phi96100.1
1 269 269
2 269 269
3 269 290
4 269 291
5 269 291
6 269 269Lo siguiente es usar la función alleles2loci() para crear el archivo tipo loci. Es importante indicar las columnas de: 1. nombres de los individuos y 2. nombres de las poblaciones. Si cada una de las columnas de los alelos tienen nombres, los nombres de la primera columna de cada genotipo se usan como nombres de los loci de salida.
En la siguiente imagen se ve cómo están estructurados los datos en el archivo maiz.csv:
Por ser individuos diploides, se coloca ploidy = 2. El argumento rownames = 1 significa que los individuos están en la columna 1, population = 2, indica que los datos de las poblaciones están en la columna 2.
alelos.maiz <- alleles2loci(data.maiz, ploidy = 2, rownames = 1, population = 2)
alelos.maiz Allelic data frame: 224 individuals
20 loci
1 additional variableSe crea un objeto que contiene una tabla de datos con los individuos, los alelos separados por “/” en una sola columna y con las poblaciones. Además de un objeto donde se almacenan todos los loci, se puede acceder a cada uno de ellos con el operador $ (por ejemplo alelos.maiz$bnlg1046).
3. Librería adegenet
Adegenet es una librería para el manejo, análisis y la simulación de datos genéticos y genómicos. Fue desarrollado por Jombart (2008) con la idea de crear una librería que conectara otros paquetes y software para el análisis de datos genéticos y de análisis multivariado.
Adegenet contiene herramientas que permiten el almacenamiento de varios tipos de datos genéticos de individuos o poblaciones. Además, pueden ser importados datos de diversos formatos generados por otras librerías.
4. Objetos genind
La mayoría de librerías que trabajan con datos genéticos, usan un tipo de objeto llamado genind. Los objetos genind son flexibles en cuanto a la información que son capaces de almacenar, se pueden incluir poblaciones, niveles jerárquicos, coordenadas geográficas, entre otros. Los objetos genind presentan los siguientes componentes:
• tab: matriz de individuos en filas y alelos en columnas
• loc.n.all: el número de alelos por cada marcador (da un rango)
• loc.fac: factor que indica a cuál columna corresponde cada marcador
• all.names: una lista del nombre de los alelos para cada locus
• ploidy: vector indicando la ploidía de cada individuo, se espera que para un mismo individuo la ploidía sea la misma, pero diferentes individuos pueden tener diferente ploidía.
• type: si los marcadores son codominantes (codom) o de presencia/ausencia (PA)
• call: comando utilizado para crear el objeto
• other: (optional): información extra en columnas (por ejemplo: fenotipos, sexo, etc)
• strata: (opcional) un data.frame que contienen la estructura jerárquica de los individuos
• pop: información sobre las poblaciones en que se encuentran los individuos
Los objetos genind se crean a partir de tablas de genotipos que deben cumplir ciertas condiciones
Los genotipos deben estar en filas , donde cada fila representa un individuo nuevo
Los loci deben estar en columnas
Se puede incluir columnas que definen la población del individuo
Los alelos deben ser caracteres separados por otro carácter como “comas”, “puntos y comas” o tabulaciones.
Para para poder trabajar con los datos en formato genind, se utiliza la función loci2genind(x, ploidy = 2, na.alleles = NA) donde se indica el objeto loci x, el nivel de ploidía y el caracter para los valores perdidos (Es importante siempre fijarse cómo están codificados los valores perdidos, en esta caso es NA, en otros puede ser 0). Con este comando que viene de la librería adegenet, convertimos el archivo a un objeto genind:
magen <- loci2genind(alelos.maiz, ploidy = 2, na.alleles = "NA")
magen /// GENIND OBJECT /////////
// 224 individuals; 20 loci; 351 alleles; size: 383.6 Kb
// Basic content
@tab: 224 x 351 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 4-37)
@loc.fac: locus factor for the 351 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: adegenet::df2genind(X = as.matrix(x[, attr(x, "locicol"), drop = FALSE]),
sep = "/", pop = pop, ploidy = ploidy)
// Optional content
@pop: population of each individual (group size range: 81-143)5. Estratificación
Cuando se trabaja con varias poblaciones, es común tener datos de sitios y agrupaciones en diferentes escalas o estratos. Es decir, se pueden tener regiones, poblaciones dentro de regiones y sub-poblaciones dentro de poblaciones. Cuando esto sucede se dice que se tiene una población estratificada jerárquicamente. El nivel de análisis dependerá de la pregunta que se desea contestar. Por ejemplo, se pueden definir las poblaciones y sub-poblaciones según cercanía geográfica, por áreas definidas políticamente o por alguna barrera natural o artificial, o por alguna característica inherente a lo que se esté estudiando (sexo, color…). Los análisis se pueden llevar a cabo en cualquiera de los estratos, lo que permite hacer comparaciones entre las distintas poblaciones, más que tener estimativas o análisis para un grupo.
Para establecer la estratificación, en adegenet, se puede utilizar la función strata (x) <- value donde x es el objeto genind y en value se indican los datos para asignar los estratos:
#Importar los estratos
estratos <- read.csv2("MaGene_Strata.csv")
strata(magen) <- estratos
nameStrata(magen) <- ~region/sitio/color #para que lo haga de manera anidada
magen /// GENIND OBJECT /////////
// 224 individuals; 20 loci; 351 alleles; size: 389.5 Kb
// Basic content
@tab: 224 x 351 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 4-37)
@loc.fac: locus factor for the 351 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: adegenet::df2genind(X = as.matrix(x[, attr(x, "locicol"), drop = FALSE]),
sep = "/", pop = pop, ploidy = ploidy)
// Optional content
@pop: population of each individual (group size range: 81-143)
@strata: a data frame with 3 columns ( region, sitio, color )6. Manejo de objetos genind
La librería adegenet incluye diversas funciones para trabajar con los objetos genind y obtener información de estos:
a- Obtener el nombre de los loci: función locNames(x)
#nombre de los loci
locNames(magen) [1] "bnlg1046" "bnlg1191" "bnlg1194" "bnlg1265" "bnlg1449" "bnlg1523"
[7] "bnlg1740" "bnlg1890" "bnlg1917" "bnlg2305" "phi015" "phi031"
[13] "phi064" "phi072" "phi078" "phi089" "phi093" "phi109188"
[19] "phi116" "phi96100" b- Obtener el número de alelos por loci: función nAll(x)
#número de alelos por loci
nAll(magen) bnlg1046 bnlg1191 bnlg1194 bnlg1265 bnlg1449 bnlg1523 bnlg1740 bnlg1890
21 24 37 28 16 34 27 32
bnlg1917 bnlg2305 phi015 phi031 phi064 phi072 phi078 phi089
26 16 10 4 16 10 10 4
phi093 phi109188 phi116 phi96100
10 9 8 9 c- Hacer selección de datos para trabajar solo con ciertos loci:
#selección para trabajar con algunos loci
m <- magen[ ,loc=c("bnlg1046","bnlg1194")]
m /// GENIND OBJECT /////////
// 224 individuals; 2 loci; 58 alleles; size: 86.8 Kb
// Basic content
@tab: 224 x 58 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 21-37)
@loc.fac: locus factor for the 58 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: .local(x = x, i = i, j = j, loc = ..1, drop = drop)
// Optional content
@pop: population of each individual (group size range: 81-143)
@strata: a data frame with 3 columns ( region, sitio, color )d- Asignar poblaciones: función setPop(x) <- formula, esta permite establecer poblaciones de trabajo. Necesario cuando se requiere realizar análisis en un estrato en específico (población, subpoblación) o una mezcla de estratos. Es decir, podemos estimar diversidad genética dentro de cada población, y luego, en un segundo análisis estimar diversidad genética en diferentes regiones. Para indicarle a la librería adegenet en cuál estrato queremos calcular los estadísticos de diversidad, se necesita definirlo con el comando setPop(x) y en la fórmula establecer el nivel jerárquico o estrato que nos interesa:
#por sitio (sitio de muestreo)
setPop(magen) <- ~sitio
magen Con la función popNames (x) se puede revisar el nombre de las poblaciones con las que se está trabajando. En el código anterior establecimos el estrato como sitio y a continuación vemos los nombres de los diferentes sitios de colecta:
popNames(magen) [1] "BORUCA" "LAUREL" "CONCEPCION" "VALENCIA" "STA_CRUZ"
[6] "LIBERIA" "CANAS" "CRUZ_LB" "CRUZ_SC" "HOJANCHA" e- Separar el objeto genind en varios objetos genind por población con seppop (x) <- formula
#comando seppop es como hacer split, separa la base de datos por poblaciones
magen2 <- seppop(magen, ~region)
magen2$BRUNCA
/// GENIND OBJECT /////////
// 81 individuals; 20 loci; 351 alleles; size: 179.9 Kb
// Basic content
@tab: 81 x 351 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 4-37)
@loc.fac: locus factor for the 351 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: .local(x = x, i = i, j = j, pop = ..1, treatOther = ..2, quiet = ..3,
drop = drop)
// Optional content
@pop: population of each individual (group size range: 81-81)
@strata: a data frame with 3 columns ( region, sitio, color )
$CHOROTEGA
/// GENIND OBJECT /////////
// 143 individuals; 20 loci; 351 alleles; size: 271 Kb
// Basic content
@tab: 143 x 351 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 4-37)
@loc.fac: locus factor for the 351 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: .local(x = x, i = i, j = j, pop = ..1, treatOther = ..2, quiet = ..3,
drop = drop)
// Optional content
@pop: population of each individual (group size range: 143-143)
@strata: a data frame with 3 columns ( region, sitio, color )f- Seleccionar datos que correspondan a un tamaño de muestra dado función selPopSize(x):
magen2 <- selPopSize(magen,pop=magen$strata$sitio,nMin=10)
magen2 /// GENIND OBJECT /////////
// 210 individuals; 20 loci; 351 alleles; size: 368.9 Kb
// Basic content
@tab: 210 x 351 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 4-37)
@loc.fac: locus factor for the 351 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: .local(x = x, i = i, j = j, drop = drop)
// Optional content
@pop: population of each individual (group size range: 12-37)
@strata: a data frame with 3 columns ( region, sitio, color )popNames(magen) [1] "BORUCA" "LAUREL" "CONCEPCION" "VALENCIA" "STA_CRUZ"
[6] "LIBERIA" "CANAS" "CRUZ_LB" "CRUZ_SC" "HOJANCHA" popNames(magen2) [1] "BORUCA" "LAUREL" "CONCEPCION" "STA_CRUZ" "LIBERIA"
[6] "CANAS" "CRUZ_LB" "CRUZ_SC" Seleccionó sólo los sitios que tengan un tamaño mínimo de 10 individuos.
7. Valores perdidos
Muchas bases de datos genéticas contienen datos perdidos, ya que hay ciertos individuos que no pueden ser genotipados correctamente. El número de datos perdidos por población es importante ya que puede sesgar de forma importante las estimativas. Se deben indicar las poblaciones con muchos valores perdidos (e.g., más de las mitad) y tomar decisiones sobre si incluirlos o no en el análisis. Para obtener los valores perdidos por marcador, se utiliza el comando info_table(x, type = "missing", percent = TRUE, plot = TRUE).
Se debe indicar el tipo de información que se desea obtener, en este caso como son valores perdidos se usa "missing". Los valores como porcentaje (default o TRUE) son específicos para cada locus, no para todos los datos en total. Si se desea que los datos se muestren como número de individuos con valores perdidos por locus, debe indicarse como percent = FALSE.
Por default lo que se obtiene es una tabla donde las poblaciones son la filas y cada marcador es una columna. Con la indicación de plot = TRUE, muestra de manera gráfica los valores por población (eje y) y por marcador (eje x).
setPop(magen) <- ~ region
perdidos <- info_table(magen, type = "missing") #por default se muestran los porcentajes y la tabla con información de datos perdidos
perdidos Locus
Population bnlg1046 bnlg1191 bnlg1194 bnlg1265 bnlg1449 bnlg1523 bnlg1740
BRUNCA 0.123 0.025 0.099 0.037 0.099 0.086 0.148
CHOROTEGA 0.175 0.105 0.133 0.091 0.203 0.105 0.140
Total 0.156 0.076 0.121 0.071 0.165 0.098 0.143
Locus
Population bnlg1890 bnlg1917 bnlg2305 phi015 phi031 phi064 phi072 phi078
BRUNCA 0.037 0.173 0.086 0.062 0.062 0.111 0.037 0.074
CHOROTEGA 0.070 0.182 0.133 0.042 0.098 0.189 0.070 0.140
Total 0.058 0.179 0.116 0.049 0.085 0.161 0.058 0.116
Locus
Population phi089 phi093 phi109188 phi116 phi96100 Mean
BRUNCA 0.012 0.062 0.025 0.074 0.012 0.072
CHOROTEGA 0.133 0.119 0.056 0.105 0.112 0.120
Total 0.089 0.098 0.045 0.094 0.076 0.103En el resultado anterior se observa que la región Chorotega tiene en promedio un 12% de valores perdidos (mean = 0.120), sobre el 7.2% de la región Brunca (Mean = 0.072). Esta cantidad de valores perdidos no es muy alta. Usualmente en microsatélites se espera que los loci tengan entre un 5 y 10% de valores perdidos. Lo importante es ver si hay algún locus que pueda ser eliminado mejorando así el promedio. Por ejemplo, el locus bnlg1917 tiene un total del 17.9% de valores perdidos. Se puede eliminar este locus y volver a correr el análisis para ver si el porcentaje de valores perdidos baja. Algunas veces es mejor eliminar un locus y mantener más individuos en el análisis. La única manera de determinar esto, es repitiendo los análisis varias veces con y sin loci que tienen muchos valores perdidos.
Para crear el gráfico donde se muestran los valores perdidos como número de individuos por marcador se puede utilizar el siguiente comando. Recuerde que si desea utilizar otro nivel de estratificación, debe utilizar el comando setPop(x) <- value.
Fig 1. Datos perdidos por locus y por población Brunca o Chorotega. El color rojo muestra un alto número de valores perdidos, mientras que el azul un número bajo de valores perdidos. Se reporta como número de individuos por marcador. Note que blng1917 tiene 40 individuos con valores perdidos.
setPop(magen) <- ~region
info_table(magen, type = "missing", percent = FALSE, plot = TRUE)
Locus
Population bnlg1046 bnlg1191 bnlg1194 bnlg1265 bnlg1449 bnlg1523 bnlg1740
BRUNCA 10.0 2.0 8.0 3.0 8.0 7.0 12.0
CHOROTEGA 25.0 15.0 19.0 13.0 29.0 15.0 20.0
Total 35.0 17.0 27.0 16.0 37.0 22.0 32.0
Locus
Population bnlg1890 bnlg1917 bnlg2305 phi015 phi031 phi064 phi072 phi078
BRUNCA 3.0 14.0 7.0 5.0 5.0 9.0 3.0 6.0
CHOROTEGA 10.0 26.0 19.0 6.0 14.0 27.0 10.0 20.0
Total 13.0 40.0 26.0 11.0 19.0 36.0 13.0 26.0
Locus
Population phi089 phi093 phi109188 phi116 phi96100 Mean
BRUNCA 1.0 5.0 2.0 6.0 1.0 5.8
CHOROTEGA 19.0 17.0 8.0 15.0 16.0 17.1
Total 20.0 22.0 10.0 21.0 17.0 23.0Fig 2. Datos perdidos por locus y por sitio. El color rojo muestra un alto número de valores perdidos, mientras que el azul un número bajo de valores perdidos. Se reporta como número de individuos por marcador. Note que blng1917 tiene 40 individuos con valores perdidos.
setPop(magen) <- ~sitio
info_table(magen, type = "missing", percent = FALSE, plot = TRUE)
Locus
Population bnlg1046 bnlg1191 bnlg1194 bnlg1265 bnlg1449 bnlg1523 bnlg1740
BORUCA 2.0 . 3.0 . 4.0 3.0 1.0
LAUREL 3.0 2.0 . . . 3.0 2.0
CONCEPCION 3.0 . 2.0 1.0 2.0 . 7.0
VALENCIA 2.0 . 3.0 2.0 2.0 1.0 2.0
STA_CRUZ 10.0 5.0 3.0 4.0 8.0 4.0 10.0
LIBERIA 3.0 3.0 4.0 3.0 10.0 3.0 1.0
CANAS 3.0 3.0 2.0 1.0 5.0 3.0 4.0
CRUZ_LB 4.0 3.0 4.0 1.0 4.0 2.0 2.0
CRUZ_SC 5.0 . 5.0 1.0 . 1.0 2.0
HOJANCHA . 1.0 1.0 3.0 2.0 2.0 1.0
Total 35.0 17.0 27.0 16.0 37.0 22.0 32.0
Locus
Population bnlg1890 bnlg1917 bnlg2305 phi015 phi031 phi064 phi072 phi078
BORUCA 2.0 5.0 4.0 2.0 . 2.0 . 2.0
LAUREL . 1.0 1.0 . . . . 1.0
CONCEPCION . 7.0 . . 2.0 3.0 2.0 1.0
VALENCIA 1.0 1.0 2.0 3.0 3.0 4.0 1.0 2.0
STA_CRUZ 5.0 6.0 6.0 2.0 3.0 7.0 3.0 7.0
LIBERIA . 6.0 3.0 3.0 6.0 9.0 3.0 4.0
CANAS 2.0 5.0 4.0 . 1.0 2.0 2.0 5.0
CRUZ_LB 1.0 3.0 3.0 1.0 2.0 6.0 1.0 2.0
CRUZ_SC 2.0 6.0 1.0 . . 1.0 . 1.0
HOJANCHA . . 2.0 . 2.0 2.0 1.0 1.0
Total 13.0 40.0 26.0 11.0 19.0 36.0 13.0 26.0
Locus
Population phi089 phi093 phi109188 phi116 phi96100 Mean
BORUCA . . . 1.0 . 1.6
LAUREL . 1.0 . . . 0.7
CONCEPCION . 2.0 2.0 3.0 1.0 1.9
VALENCIA 1.0 2.0 . 2.0 . 1.7
STA_CRUZ 5.0 5.0 1.0 4.0 6.0 5.2
LIBERIA 5.0 4.0 4.0 3.0 3.0 4.0
CANAS 1.0 3.0 2.0 1.0 4.0 2.6
CRUZ_LB 5.0 2.0 . 4.0 2.0 2.6
CRUZ_SC 2.0 1.0 1.0 . . 1.4
HOJANCHA 1.0 2.0 . 3.0 1.0 1.2
Total 20.0 22.0 10.0 21.0 17.0 23.0Existen diversas opciones para almacenar los datos de valores perdidos o de diversidad por locus. Una de éstas es el comando write.csv(x, file = "nombre.csv"), donde “x” es el objeto de los valores perdidos y en “file =” se indica el nombre del archivo de salida, debe ir entre comillas y con la extensión .csv. Deben verificar en la carpeta correspondiente si se creó el archivo.
write.csv(perdidos, file= "Valores_perdidos.csv") 8. Importar datos desde un VCF
Un archivo VCF puede dividirse en tres secciones: un encabezado VCF, una región fija y una región gt (genotipos). La región meta del VCF se encuentra en la parte superior del archivo y contiene metadatos que describen el cuerpo del archivo. Cada línea meta del VCF comienza con ‘##’. La información en la región meta define las abreviaturas utilizadas en otras partes del archivo y puede documentar el software utilizado para crear el archivo, así como los parámetros utilizados por este software.
Fig 3. Estructura de un archivo VCF
Para importar archivos desde un VCF, se debe instalar y cargar el paquete vcfR (Knaus & Grünwald, 2017)
#install.packages("vcfR")
library(vcfR)
pajaros <- read.vcfR("Pajaros.vcf")Scanning file to determine attributes.
File attributes:
meta lines: 10
header_line: 11
variant count: 1925
column count: 49
Meta line 10 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
Character matrix gt rows: 1925
Character matrix gt cols: 49
skip: 0
nrows: 1925
row_num: 0
Processed variant 1000
Processed variant: 1925
All variants processedpajaros***** Object of Class vcfR *****
40 samples
594 CHROMs
1,925 variants
Object size: 1.6 Mb
0 percent missing data
***** ***** *****gen.pajaros <- vcfR2genind(pajaros) # gen.pajaros es un objeto genind.
gen.pajaros/// GENIND OBJECT /////////
// 40 individuals; 1,925 loci; 3,850 alleles; size: 1.8 Mb
// Basic content
@tab: 40 x 3850 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 2-2)
@loc.fac: locus factor for the 3850 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: adegenet::df2genind(X = t(x), sep = sep)
// Optional content
- empty -# ahora crear estratos para el objeto gen.pajaros
# importamos la tabla de datos que tiene la definicion de poblacion
# como esta separado por tabulaciones, usamos comando read.delim()
estratos.pajaros <- read.delim("Sporo_pop.txt")
head(estratos.pajaros) ID poblacion
1 Sc-coDO-727 Sarapiqui
2 Sc-coDO-743 Sarapiqui
3 Sc-coDO-744 Sarapiqui
4 Sc-coUCR042 Sarapiqui
5 Sc-coUCR077 Sarapiqui
6 Sc-coUM-014 Sarapiqui## hay que asignar la jerarquia de estratos al objeto genind.
## para que los ananlisis se hagan sobre la estructura correcta
## esto se hace con el comando strata() <- marco de datos con los estratos
strata(gen.pajaros) <- estratos.pajaros
gen.pajaros/// GENIND OBJECT /////////
// 40 individuals; 1,925 loci; 3,850 alleles; size: 1.8 Mb
// Basic content
@tab: 40 x 3850 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 2-2)
@loc.fac: locus factor for the 3850 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: adegenet::df2genind(X = t(x), sep = sep)
// Optional content
@strata: a data frame with 2 columns ( ID, poblacion )## ahora definimos sobre cual de los estratos se realizaran los analisis
## esto se hace con el comando setPop() que define cual jerarquia de strata
## se utilizara. Del comando se coloca una formula que indica la
## estructura jerarquica. En este caso solo hay *una sola* estrato, poblacion
setPop(gen.pajaros) <- ~poblacion
## para corroborar, vemos la poblacion asignada a cda individuo
gen.pajaros@pop [1] Sarapiqui Sarapiqui Sarapiqui Sarapiqui Sarapiqui Sarapiqui Sarapiqui
[8] Sarapiqui Sarapiqui Sarapiqui Cahuita Cahuita Cahuita Cahuita
[15] Cahuita Cahuita Cahuita Cahuita Cahuita Cahuita Golfito
[22] Golfito Golfito Golfito Golfito Golfito Golfito Golfito
[29] Golfito Golfito Quepos Quepos Quepos Quepos Quepos
[36] Quepos Quepos Quepos Quepos Quepos
Levels: Sarapiqui Cahuita Golfito Quepos9. Exportar datos
La función genind2df() permite la conversión de objetos genind a una tabla de datos para trabajarla en otras bases de datos o para ser exportada a un archivo de texto.
maiz <- genind2df(magen2, usepop = T, oneColPerAll = T) #para tener un alelo por columna, si se quisieran unirlos por algún se utiliza el argumento sep = "/" Referencias consultadas
Bedoya, C. (2012). Estudios de diversidad genética en poblaciones de maíz (Zea mays L.) evaluadas con microsatélites. Tesis doctoral. Universidad de las Islas Baleares, España. 223p.
Fuchs, E. & Barrantes, G. (2015). El lenguaje estadístico R aplicado a las ciencias biológicas. Editorial UCR, San José, Costa Rica. 197p.
Jiménez, R. (2014). Caracterización de las razas criollas e indígenas de maíz colombiano por medio de marcadores moleculares SSR. Tesis para optar por el grado de Magister en Ciencias Biológicas. Universidad Nacional de Colombia Facultad de Ciencias Agropecuarias. 70p.
Grünwald, N., Kamvar, Z. & Everhart, S. (2015). Population genetics in R. Corvallis, Oregon, USA. Recuperado de http://grunwaldlab.github.io/Population_Genetics_in_R/Data_Preparation.html (Consulta 02 febrero 2017)
Paradis E. 2010. pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics 26: 419-420.
Grunwald Lab. Reading VCF Files. Recuperado de https://grunwaldlab.github.io/Population_Genetics_in_R/reading_vcf.html (Consulta 20 Marzo 2024)
R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.